%%   导航信号跟踪演示模型
%   Copyright 2022, NUDT
%
%   国防科技大学电子科学学院
%   导航与时空技术工程研究中心
%
%   说明：
%       多径干扰信号误差分析
%
%   输入：
%       无
%   输出：
%       无
%
%   原始作者：刘吉洋
%   建立日期：2023年11月17日
%
%   更新历史:
%
%% 简单多路径仿真
alapha=0.5;
theat=0:0.01:2*pi;
y_end=sin(theat);
error1=atan(alapha*sin(theat)./(1+alapha*cos(theat)));
plot(theat,error1)
xlabel('多径相位');
ylabel('载波环相位误差');
title ('载波相位误差');
grid on
%% 平均效果仿真
c=3e+8;
a=pi/3;
b=pi/6:0.01:5*pi/6;
delta_tau1=2*5*sin(a+b)/c; %多径初始时延
fc=1575.42e6;%载波频率
delta_fai_1=delta_tau1*2*pi*fc*1e-6;%多径相位时延变化
alapha1=0.5;%衰落系数
fai_c=atan(alapha1*sin(delta_fai_1)./(1+alapha1*cos(delta_fai_1)));%载波跟踪环路相位误差
fai=pi;%相位
Tc=1;%C/A码片宽度/s
d=Tc/2;%早迟间隔
tau_L=(1+alapha1*cos(fai))/2*d;
tau_H=(alapha1*d*cos(fai))/2+Tc-d/2;
error1=(alapha1.*delta_tau1.*cos(fai)./(1+alapha1.*cos(fai))).*(delta_tau1<=tau_L)+alapha1.*d/2*cos(fai).*(delta_tau1>tau_L & delta_tau1<tau_H )+alapha1.*cos(fai)./(2-alapha1.*cos(fai)).*(d/2+Tc-delta_tau1).*(delta_tau1>=tau_H & delta_tau1<Tc+d/2);
hold on
plot(delta_tau1,error1,'-');
xlabel('多径时延 / 码片');
ylabel('码环误差/码片');
title ('DLL误差包络');
grid on
figure
%plot(delta_tau1,fai_c);
xlabel('多径时延 / 码片');
ylabel('载波环相位误差');
title ('DLL误差包络');
grid on
%% 多条多径误差验证
clear
M=10;
for ll=1:M
    tau_number = rand(1,ll);
    tau = tau_number * 1;%时延
    alapha_number = rand(1,ll);
    alapha = alapha_number * 0.1;%反射系数
    error_theory=0;
    for ii =1:length(tau)
        delta_tau1=tau(ii); %多径初始时延
        fc=1575.42e6;%载波频率
        delta_fai_1=delta_tau1*2*pi*fc*1e-6;%多径相位时延变化
        alapha1=alapha(ii);%衰落系数
        %fai_c=atan(alapha1*sin(delta_fai_1)./(1+alapha1*cos(delta_fai_1)));%载波跟踪环路相位误差
        fai=pi;%相位
        Tc=1;%C/A码片宽度/s
        d=0.5;%早迟间隔
        tau_L=(1+alapha1*cos(fai))/2*d;
        tau_H=(alapha1*d*cos(fai))/2+Tc-d/2;
        error1=(alapha1.*delta_tau1.*cos(fai)./(1+alapha1.*cos(fai))).*(delta_tau1<=tau_L)+alapha1.*d/2*cos(fai).*(delta_tau1>tau_L & delta_tau1<tau_H )+alapha1.*cos(fai)./(2-alapha1.*cos(fai)).*(d/2+Tc-delta_tau1).*(delta_tau1>=tau_H & delta_tau1<Tc+d/2);
        error_theory=error_theory+error1;
    end
    error_theory_M(ll)=error_theory;

    x=-2:0.001:2;
    y=(x+1).*(x<=0 & x>-1)-(x-1).*(x<1 & x>0);
    y_end=0;
    for ii =1:length(tau)
        tau1=tau(ii)*1000;
        alapha1=alapha(ii);
        y1=zeros(1,length(y));
        y1(tau1+1:end)=y(1:end-tau1);
        y1(1:tau1)=y(end-tau1+1:end);
        y_end=y_end-alapha1*y1;
    end
    y_end=y_end+y;

    m=0;
    f=find(y_end==max(y_end));
    d=500;%相关距离
    index=10;
    while abs(index+m-d)>2
        m=m+1;
        i=1;
        min=10;
        while i<1000
            valu=y_end(f-m);
            cha=abs(valu-y_end(f+i));
            if cha<min
                min=cha;
                index=i;
                value1=y_end(f-m);
                value2=y_end(f+i);
            end
            i=i+1;
        end
    end
    error_sim=(index-m)/2/1000;
    error_sim_M(ll)=error_sim;
end
figure
plot(x,y_end,'-')
set(gca,'XLim',[-2 2]);%X轴的数据显示范围
set(gca,'YLim',[0 1.5]);%X轴的数据显示范围
figure
plot(error_sim_M,'o')
hold on
plot(error_theory_M,'-')
%% 抖动码环误差仿真
clear;
delta_tau0=0:0.08:2; %多径初始时延
for ll=1:length(delta_tau0)
    fc=1575.42e6;%载波频率
    L=0.001;%时延变化范围（码片）
    delta_t=-L:L/400:L;%抖动时延
    delta_tau1=delta_tau0(ll)+delta_t;
    alapha1=0.5;%衰落系数
    delta_fai_1=delta_tau1*2*pi*fc*1e-6;%多径相位时延变化
    delta_fai=mod(delta_fai_1,2*pi);
    fai_c=atan(alapha1*sin(delta_fai_1)./(1+alapha1*cos(delta_fai_1)));%载波跟踪环路相位
    fai=delta_fai_1;%相位
    Tc=1;%C/A码片宽度
    d=Tc/8;%早迟间隔
    tau_L=(1+alapha1*cos(fai))/2*d;
    tau_H=(alapha1*d*cos(fai))/2+Tc-d/2;
    error1=(alapha1.*delta_tau1.*cos(fai)./(1+alapha1.*cos(fai))).*(delta_tau1<=tau_L)+alapha1.*d/2.*cos(fai).*(delta_tau1>tau_L & delta_tau1<tau_H )+alapha1.*cos(fai)./(2-alapha1.*cos(fai)).*(d/2+Tc-delta_tau1).*(delta_tau1>tau_H & delta_tau1<Tc+d/2);
    error_Sim(ll)=mean(error1);
    k=-0.9716;
end
%理论结果
alapha1=0.5;%衰落系数
%fai_c=atan(alapha1*sin(delta_fai_1)./(1+alapha1*cos(delta_fai_1)));%载波跟踪环路相位误差
k2=0.2227;
delta_tau01=0:0.005:2; %多径初始时延
tau_L1=(1-alapha1)/2*d;
tau_L2=(1+alapha1)/2*d;
tau_H1=-(alapha1*d)/2+Tc-d/2;
tau_H2=(alapha1*d)/2+Tc-d/2;
error_Ana=(delta_tau01*k/(2*pi)).*(delta_tau01<=tau_L1)+...
(-(1-alapha1)*delta_tau01*k/(4*pi*alapha1)+(1-alapha1^2)*d*k/(8*pi*alapha1)).*(delta_tau01<=tau_L2 & delta_tau01>tau_L1)+...
0*(delta_tau01>tau_L2 & delta_tau01<tau_H1)-...
k2*(alapha1-2)./(4*pi*alapha1).*(delta_tau01-Tc+(alapha1+1)*d/2).*(delta_tau01>tau_H1 & delta_tau01<=tau_H2)+...
k2./(2*pi).*(d/2+Tc-delta_tau01).*(delta_tau01>tau_H2 & delta_tau01<Tc+d/2);
hold on
plot(delta_tau0,error_Sim,'>')
hold on
plot(delta_tau01,error_Ana,'-')

%%
hold on
grid on
xlim([0,1.6])
ylim([-0.05,0.04])
xlabel('多径时延 / 码片');
ylabel('码相位误差 / 码片');
%% 抖动随抖动幅度变化
clear;
delta_tau0=0.042; %多径初始时延（码片）
l=0.0001:0.00008:0.005;
kesai=0;
fc=1575.42e6;%载波频率
r=1/(fc*1e-6)*2.93/3;
LL=abs(l.*cos(kesai));%时延变化范围（码片）
for ll=1:length(LL)
    L=LL(ll);
    delta_t=0:L/100:L;%抖动时延
    delta_tau1=delta_tau0+delta_t;
    alapha1=0.5;%衰落系数
    fai_0=3;%载波初相
    delta_fai_1=delta_tau1*2*pi*fc*1e-6;%多径相位时延变化
    fai_c=fai_0+1+atan(alapha1*sin(delta_fai_1+fai_0)./(1+alapha1*cos(delta_fai_1+fai_0)));%载波跟踪环路相位
    fai=delta_fai_1;%相位
    Tc=1;%C/A码片宽度
    d=Tc/2;%早迟间隔
    tau_L=(1+alapha1*cos(fai))/2*d;
    tau_H=(alapha1*d*cos(fai))/2+Tc-d/2;
    error1=(alapha1.*delta_tau1.*cos(fai)./(1+alapha1.*cos(fai))).*(delta_tau1<=tau_L)+alapha1.*d/2.*cos(fai).*(delta_tau1>tau_L & delta_tau1<tau_H )+alapha1.*cos(fai)./(2-alapha1.*cos(fai)).*(d/2+Tc-delta_tau1).*(delta_tau1>tau_H & delta_tau1<Tc+d/2);
    error_sim(ll)=abs(mean(error1));
    % 理论结果
    a=alapha1;
    x=fai;
    x_begin=mod(x(1)-pi,2*pi)+pi;
    x_end=mod(x(end)-pi,2*pi)+pi;
    aa=x(1)-x_begin;
    bb=-x(end)+x_end;
    d=(aa+bb)/2/pi;
    y_1=log((a+cos(x_begin)+sin(x_begin)*sqrt(a^2-1))./(a*cos(x_begin)+1))./(sqrt(a^2-1));
    y_end=log((a+cos(x(end))+sin(x_end)*sqrt(a^2-1))./(a*cos(x_end)+1))./(sqrt(a^2-1));
    error_anl(ll)=delta_tau0*(-y_end+y_1+x_end-x_begin+d*0.9716)./(x(end)-x(1));
end
lll=0.0001:0.000008:0.005;
for ll=1:length(lll)
    L=lll(ll);
    delta_t=0:L/100:L;%抖动时延
    delta_tau1=delta_tau0+delta_t;
    alapha1=0.5;%衰落系数
    fai_0=3;%载波初相
    delta_fai_1=delta_tau1*2*pi*fc*1e-6;%多径相位时延变化
    fai_c=fai_0+1+atan(alapha1*sin(delta_fai_1+fai_0)./(1+alapha1*cos(delta_fai_1+fai_0)));%载波跟踪环路相位
    fai=delta_fai_1;%相位
    % 理论结果
    a=alapha1;
    x=fai;
    x_begin=mod(x(1)-pi,2*pi)+pi;
    x_end=mod(x(end)-pi,2*pi)+pi;
    aa=x(1)-x_begin;
    bb=-x(end)+x_end;
    d=(aa+bb)/2/pi;
    y_1=log((a+cos(x_begin)+sin(x_begin)*sqrt(a^2-1))./(a*cos(x_begin)+1))./(sqrt(a^2-1));
    y_end=log((a+cos(x(end))+sin(x_end)*sqrt(a^2-1))./(a*cos(x_end)+1))./(sqrt(a^2-1));
    error_anl(ll)=delta_tau0*(-y_end+y_1+x_end-x_begin+d*0.9716)./(x(end)-x(1));
end
hold on
plot(l,error_sim,'>')
plot(lll,abs(error_anl),'-')

%%
hold on
xlim([0.0003,0.003])
ylim([0,0.02])
xlabel('抖动幅度 / 码片');
ylabel('误差 / 码片');
title ('DLL误差');
grid on
%% 理论最佳
hold on
y_end=0:0.02:0.02;
r=zeros(length(y_end));
r=r+1/1575*3/2.93;
plot(4*r,y_end,'-')
%%
xlabel('多径时延 / 码片');
ylabel('码环误差/码片');
ylim([-0.25,0.25]);
xlim([0.07,0.15]);
%%
hold on
plot(delta_tau0,error,'-->');
xlabel('多径时延 / 码片');
ylabel('码环误差/码片');
axis([0 2 -2.5 2.5]);
title ('DLL误差包络');
grid on
%% 抖动最佳条件下多条多径情况
clear;
m=1:1:30; %多径条数
fc=1575.42e6; %载波频率
sta=pi/6;%卫星仰角
la1=pi/3;%反射面倾角
c=3e+8;%光速
r=10*c/(4*fc*sin(la1+sta));%抖动幅度(cm)
L=4.4*r*sin(la1+sta)/c*1e6;%时延变化范围
delta_tau0=0.01; %多径初始时延（码片）
N=100;
delta_t=0:L/N:L;%直线抖动时延
delta_tau1=delta_tau0+delta_t;
alapha1=0.5;%衰落系数
delta_fai_1=delta_tau1*2*pi*fc*1e-6;%多径相对相位变化
fai=delta_fai_1;%相位
Tc=1;%C/A码片宽度
d=Tc/8;%早迟间隔
tau_L=(1+alapha1*cos(fai))/2*d;
tau_H=(alapha1*d*cos(fai))/2+Tc-d/2;
error1=(alapha1.*delta_tau1.*cos(fai)./(1+alapha1.*cos(fai))).*(delta_tau1<=tau_L)+alapha1.*d/2.*cos(fai).*(delta_tau1>tau_L & delta_tau1<tau_H )+alapha1.*cos(fai)./(2-alapha1.*cos(fai)).*(d/2+Tc-delta_tau1).*(delta_tau1>tau_H & delta_tau1<Tc+d/2);
error1_main=mean(error1)
Q=100;%蒙特卡洛仿真次数
for pp=1:Q
    for ll=2:length(m)
        M=m(ll);
        error2=0;
        for dd=1:M
            la_number = rand;
            la = la_number * pi/2 + pi/4;%反射面倾角
            delta_tau0_number = rand;
            delta_tau0 = delta_tau0_number * 2 + 0;%多径初始时延（码片）
            p=abs(la1-la);% 抖动角度
            L=4.4*r*cos(p)*sin(la1+sta)/c*1e6; %相对时延变化范围
            delta_t=0:L/N:L;%直线抖动时延
            delta_tau1=delta_tau0+delta_t;
            alapha1_number = rand;
            alapha1 = delta_tau0_number * 0.1 + 0;%反射系数
            delta_fai_1=delta_tau1*2*pi*fc*1e-6;%多径相对相位变化
            fai=delta_fai_1;%相位
            Tc=1;%C/A码片宽度
            d=Tc/8;%早迟间隔
            tau_L=(1+alapha1*cos(fai))/2*d;
            tau_H=(alapha1*d*cos(fai))/2+Tc-d/2;
            error1=(alapha1.*delta_tau1.*cos(fai)./(1+alapha1.*cos(fai))).*(delta_tau1<=tau_L)+alapha1.*d/2.*cos(fai).*(delta_tau1>tau_L & delta_tau1<tau_H )+alapha1.*cos(fai)./(2-alapha1.*cos(fai)).*(d/2+Tc-delta_tau1).*(delta_tau1>tau_H & delta_tau1<Tc+d/2);
            error2(dd)=mean(error1);
        end
        error(pp,ll)=sum(error2);
    end
    error_middle=error+error1_main;
end
error_mean=mean(error_middle);
error_std=std(error_middle);
hold on
x = m;                     
y_end = abs(error_mean); %误差均值​
xconf = [x x(end:-1:1)];%构造正反向的x值，作为置信区间的横坐标值
yconf = [y_end+error_std y_end(end:-1:1)-error_std(end:-1:1)];%构造y方向的上下范围，作为置信区间的纵坐标值​
hold on
p = fill(xconf,yconf,'red');%定义填充区间
p.FaceColor = [1 0.8 0.8];%定义区间的填充颜色      
p.EdgeColor = 'none';%定义区间边界的填充颜色，此处不设置​
hold on
plot(x,y_end,'ro')%绘制折线
xlabel('多径数量');
ylabel('码环误差/码片');
grid on
%% 非最佳抖动最佳下多条多径情况
clear;
m=1:1:30; %多径条数
fc=1575.42e6; %载波频率
sta=pi/6;%卫星仰角
la1=pi/3;%反射面倾角
c=3e+8;%光速
r=1.5*c/(4*fc*sin(la1+sta));%抖动幅度(cm)
p_number= rand;
ppp = p_number * pi/2 + pi/4;%反射面倾角
L=4.4*cos(ppp)*r*sin(la1+sta)/c*1e6;%时延变化范围
delta_tau0=0.01; %多径初始时延（码片）
N=100;
delta_t=0:L/N:L;%直线抖动时延
delta_tau1=delta_tau0+delta_t;
alapha1=0.5;%衰落系数
delta_fai_1=delta_tau1*2*pi*fc*1e-6;%多径相对相位变化
fai=delta_fai_1;%相位
Tc=1;%C/A码片宽度
d=Tc/8;%早迟间隔
tau_L=(1+alapha1*cos(fai))/2*d;
tau_H=(alapha1*d*cos(fai))/2+Tc-d/2;
error1=(alapha1.*delta_tau1.*cos(fai)./(1+alapha1.*cos(fai))).*(delta_tau1<=tau_L)+alapha1.*d/2.*cos(fai).*(delta_tau1>tau_L & delta_tau1<tau_H )+alapha1.*cos(fai)./(2-alapha1.*cos(fai)).*(d/2+Tc-delta_tau1).*(delta_tau1>tau_H & delta_tau1<Tc+d/2);
error1_main=mean(error1)
Q=100;%蒙特卡洛仿真次数
for pp=1:Q
    for ll=2:length(m)
        M=m(ll);
        error2=0;
        for dd=1:M
            la_number = rand;
            la = la_number * pi/2 + pi/4;%反射面倾角
            delta_tau0_number = rand;
            delta_tau0 = delta_tau0_number * 2 + 0;%多径初始时延（码片）
            p=abs(la1-la);% 抖动角度
            L=4.4*r*cos(p)*sin(la1+sta)/c*1e6; %相对时延变化范围
            delta_t=0:L/N:L;%直线抖动时延
            delta_tau1=delta_tau0+delta_t;
            alapha1_number = rand;
            alapha1 = delta_tau0_number * 0.1 + 0;%反射系数
            delta_fai_1=delta_tau1*2*pi*fc*1e-6;%多径相对相位变化
            fai=delta_fai_1;%相位
            Tc=1;%C/A码片宽度
            d=Tc/8;%早迟间隔
            tau_L=(1+alapha1*cos(fai))/2*d;
            tau_H=(alapha1*d*cos(fai))/2+Tc-d/2;
            error1=(alapha1.*delta_tau1.*cos(fai)./(1+alapha1.*cos(fai))).*(delta_tau1<=tau_L)+alapha1.*d/2.*cos(fai).*(delta_tau1>tau_L & delta_tau1<tau_H )+alapha1.*cos(fai)./(2-alapha1.*cos(fai)).*(d/2+Tc-delta_tau1).*(delta_tau1>tau_H & delta_tau1<Tc+d/2);
            error2(dd)=mean(error1);
        end
        error(pp,ll)=sum(error2);
    end
end
error_middle=error+error1_main;
error_mean=mean(error_middle);
error_std=std(error_middle);
hold on
x = m;                     
y_end = abs(error_mean); %误差均值​
xconf = [x x(end:-1:1)];%构造正反向的x值，作为置信区间的横坐标值
yconf = [y_end+error_std y_end(end:-1:1)-error_std(end:-1:1)];%构造y方向的上下范围，作为置信区间的纵坐标值​
hold on
p = fill(xconf,yconf,'blue');%定义填充区间
p.FaceColor = [0.8 0.8 1];%定义区间的填充颜色      
p.EdgeColor = 'none';%定义区间边界的填充颜色，此处不设置​
hold on
plot(x,y_end,'ro')%绘制折线
xlabel('多径数量');
ylabel('码环误差/码片');
grid on
%% 天线静止条件下多条多径情况
clear;
m=1:1:30; %多径条数
fc=1575.42e6; %载波频率
sta=pi/6;%卫星仰角
la1=pi/3;%反射面倾角
c=3e+8;%光速
delta_tau0=0.015; %多径初始时延（码片）
delta_tau1=delta_tau0;
alapha1=0.5;%衰落系数
delta_fai_1=delta_tau1*2*pi*fc*1e-6;%多径相对相位变化
fai=delta_fai_1;%相位
Tc=1;%C/A码片宽度
d=Tc/3;%早迟间隔
tau_L=(1+alapha1*cos(fai))/2*d;
tau_H=(alapha1*d*cos(fai))/2+Tc-d/2;
error1=(alapha1.*delta_tau1.*cos(fai)./(1+alapha1.*cos(fai))).*(delta_tau1<=tau_L)+alapha1.*d/2.*cos(fai).*(delta_tau1>tau_L & delta_tau1<tau_H )+alapha1.*cos(fai)./(2-alapha1.*cos(fai)).*(d/2+Tc-delta_tau1).*(delta_tau1>tau_H & delta_tau1<Tc+d/2);
error1_main=mean(error1)
Q=100;%蒙特卡洛仿真次数
for q=1:Q
    for ll=2:length(m)
        M=m(ll);
        error2=0;
        for dd=1:M
            la_number = rand;
            la = la_number * pi/2 + pi/4;%反射面倾角
            delta_tau0_number = rand;
            delta_tau0 = delta_tau0_number * 2 + 0;%多径初始时延（码片）
            delta_tau1=delta_tau0;
            alapha1_number = rand;
            alapha1 = delta_tau0_number * 0.1 + 0;%反射系数
            delta_fai_1=delta_tau1*2*pi*fc*1e-6;%多径相对相位变化
            fai=delta_fai_1;%相位
            Tc=1;%C/A码片宽度
            d=Tc/8;%早迟间隔
            tau_L=(1+alapha1*cos(fai))/2*d;
            tau_H=(alapha1*d*cos(fai))/2+Tc-d/2;
            error1=(alapha1.*delta_tau1.*cos(fai)./(1+alapha1.*cos(fai))).*(delta_tau1<=tau_L)+alapha1.*d/2.*cos(fai).*(delta_tau1>tau_L & delta_tau1<tau_H )+alapha1.*cos(fai)./(2-alapha1.*cos(fai)).*(d/2+Tc-delta_tau1).*(delta_tau1>tau_H & delta_tau1<Tc+d/2);
            error2(dd)=mean(error1);
        end
        error(q,ll)=sum(error2);
    end
    error_middle=error+error1_main;
end
error_mean=mean(error_middle);
error_std=std(error_middle);
hold on
x = m;                     
y_end = abs(error_mean); %误差均值​
xconf = [x x(end:-1:1)];%构造正反向的x值，作为置信区间的横坐标值
yconf = [y_end+error_std y_end(end:-1:1)-error_std(end:-1:1)];%构造y方向的上下范围，作为置信区间的纵坐标值​
hold on
p = fill(xconf,yconf,'red');%定义填充区间
p.FaceColor = [1 0.8 0.8];%定义区间的填充颜色      
p.EdgeColor = 'none';%定义区间边界的填充颜色，此处不设置​
hold on
plot(x,y_end,'ro')%绘制折线
xlabel('多径数量');
ylabel('码环误差/码片');
grid on
%% 抖动误差随幅度变化
clear;
delta_tau0=0.13; %多径初始时延（码片）
[l,kesai]=meshgrid(0.0001:0.00001:0.001,0:0.01:2*pi);
LL=abs(l.*cos(kesai));%时延变化范围（码片）
[a,b]=size(LL);
for ll=1:a
    for hh=1:b %仿真
    fc=1575.42e6;%载波频率
    L=LL(ll,hh);
    delta_t=0:L/100:L;%抖动时延
    delta_tau1=delta_tau0+delta_t;
    alapha1=0.5;%衰落系数
    fai_0=3;%载波初相
    delta_fai_1=delta_tau1*2*pi*fc*1e-6;%多径相位时延变化
    fai_c=fai_0+1+atan(alapha1*sin(delta_fai_1+fai_0)./(1+alapha1*cos(delta_fai_1+fai_0)));%载波跟踪环路相位
    fai=delta_fai_1+fai_c;%相位
    Tc=1;%C/A码片宽度
    d=Tc/2;%早迟间隔
    tau_L=(1+alapha1*cos(fai))/2*d;
    tau_H=(alapha1*d*cos(fai))/2+Tc-d/2;
    error1=(alapha1.*delta_tau1.*cos(fai)./(1+alapha1.*cos(fai))).*(delta_tau1<=tau_L)+alapha1.*d/2.*cos(fai).*(delta_tau1>tau_L & delta_tau1<tau_H )+alapha1.*cos(fai)./(2-alapha1.*cos(fai)).*(d/2+Tc-delta_tau1).*(delta_tau1>tau_H & delta_tau1<Tc+d/2);
    error2(ll,hh)=abs(mean(error1));
    end
end
% 理论结果
hold on
plot(l,error2)
meshc(l,kesai/pi,error2)
xlabel('抖动幅度 / 码片');
ylabel('抖动角度 / pi');
title ('DLL误差');
grid on
%%
hold on
ylim([0,2])
%%
[x,y_end]=meshgrid(-4:0.1:4,-3:0.1:3);
z=(x.^2-2*x).*exp(-x.^2-y_end.^2-x.*y_end);
meshc(x,y_end,z)
xlabel('抖动幅度 / 码片');
ylabel('码环误差/码片');
%% 未抖动误差随时间变化
delta_tau1=0.5; %多径初始时延
t=0:1:1000;%ms
fc=1575.42e6;%载波频率
fai_0=3;%载波初相
delta_fai_1=delta_tau1*2*pi*fc*1e-6;%多径相位时延变化
fai_c=atan(alapha1*sin(delta_fai_1+fai_0)./(1+alapha1*cos(delta_fai_1+fai_0)));%载波跟踪环路相位
fai=delta_fai_1;%相位
alapha1=0.5;%衰落系数
Tc=1;%C/A码片宽度/s
d=Tc;%早迟间隔
tau_L=(1+alapha1*cos(fai))/2*d;
tau_H=(alapha1*d*cos(fai))/2+Tc-d/2;
error1=(alapha1.*delta_tau1.*cos(fai)./(1+alapha1.*cos(fai))).*(delta_tau1<=tau_L)+alapha1.*d./2*cos(fai).*(delta_tau1>tau_L & delta_tau1<tau_H )+alapha1.*cos(fai)./(2-alapha1.*cos(fai)).*(d/2+Tc-delta_tau1).*(delta_tau1>tau_H & delta_tau1<Tc+d/2);
error=abs(error1)+0.02*randn(1,length(t));
hold on
plot(t,error,'-');
xlabel('多径时延 / 码片');
ylabel('码环误差/码片');
title ('DLL误差');
grid on
%% 光线不平行入射时
H=20000e3;%同步轨道卫星高度/m
L=sqrt(3)*H/2+1000;%卫星投影点到天线距离/m
h=3;%天线高度/m
tau=sqrt((H+h).^2+L^2)-sqrt((H-h).^2+L^2);%光程误差/m
fai=2*pi*tau/(3e8);
%%
clear
Tc=1;%C/A码片宽度/s
tau=-2:0.001:2;%相关时延
R_tau=(-1/Tc*abs(tau)+1).*(tau<Tc & -Tc<tau);%直达信号自相关函数
R_tau1=(-1/Tc*abs(tau-0.5)+1).*(tau-0.5<Tc & -Tc<tau-0.5);%多径信号自相关函数
R_tau2=(-1/Tc*abs(tau-0.6)+1).*(tau-0.6<Tc & -Tc<tau-0.6);%多径信号自相关函数
figure
subplot(2,1,1)
plot(tau,R_tau)
hold on
plot(tau,0.5*R_tau1)
plot(tau,0.5*R_tau1+R_tau)
xlabel('Time delay (chips)');
ylabel('Correlation');
title('Relative phase=0°')
grid on
subplot(2,1,2)
plot(tau,R_tau)
hold on
plot(tau,R_tau)
plot(tau,-0.5*R_tau1)
plot(tau,-0.5*R_tau1+R_tau)
xlabel('Time delay (chips)');
ylabel('Correlation');
title('Relative phase=180°')
grid on
%% 随时间变化特性
delta_tau1=0.5; %多径初始时延
t=0:1:1000;%ms
fc=1575.42e6;%载波频率
fai_0=pi;%载波初相
delta_fai_1=delta_tau1*2*pi*fc*1e-6;%多径相位时延变化
fai_c=atan(alapha1*sin(delta_fai_1+fai_0)./(1+alapha1*cos(delta_fai_1+fai_0)));%载波跟踪环路相位
fai=fai_0+delta_fai_1-fai_c;%相位
alapha1=0.5;%衰落系数
Tc=1;%C/A码片宽度/s
d=Tc;%早迟间隔
tau_L=(1+alapha1*cos(fai))/2*d;
tau_H=(alapha1*d*cos(fai))/2+Tc-d/2;
error1=(alapha1.*delta_tau1.*cos(fai)./(1+alapha1.*cos(fai))).*(delta_tau1<=tau_L)+alapha1.*d./2*cos(fai).*(delta_tau1>tau_L & delta_tau1<tau_H )+alapha1.*cos(fai)./(2-alapha1.*cos(fai)).*(d/2+Tc-delta_tau1).*(delta_tau1>tau_H & delta_tau1<Tc+d/2);
error=error1+0.02*randn(1,length(t));
figure
plot(t,error,'-');
xlabel('时间 / ms');
ylabel('码环误差/码片');
title ('DLL误差');
grid on
%% 自相关
tau1=6; %时延/码片
tau2=6; %时延/码片
alpha1=0.5;%衰落系数
alpha2=0.3;%衰落系数
prncodeBase0=genCode;%直达信号
prncodeBase1=zeros(1,length(prncodeBase0));
prncodeBase1(1+tau1:end)=prncodeBase0(1:end-tau1);
prncodeBase1(1:tau1)=prncodeBase0(end-tau1+1:end);

prncodeBase2=zeros(1,length(prncodeBase0));
prncodeBase2(1+tau2:end)=prncodeBase0(1:end-tau2);
prncodeBase2(1:tau2)=prncodeBase0(end-tau2+1:end);

prncodeBase_r=prncodeBase0+alpha1*prncodeBase1+alpha2*prncodeBase2;
[c,lags] = xcorr(prncodeBase_r,prncodeBase0,'normalized');
lag=lags;
i=1;
m=0;
f=find(c==max(c));
d=9;
min=10;
while index+m~=d
    m=m+1;
    i=1;
    min=10;
    while i<8
        valu=c(f-m);
        cha=abs(valu-c(f+i));
        if cha<min
            min=cha;
            index=i;
            value1=c(f-m);
            value2=c(f+i);
        end
        i=i+1;
    end
end
error=(index+m)/2
plot(lag/8,c,'o')
set(gca,'XLim',[-2 2]);%X轴的数据显示范围
set(gca,'YLim',[0 1.5]);%X轴的数据显示范围
%% 采样扩频码
t=(1:genSampleCount) / sampleFreq;
genCode1=prncodeBase(1:1000);
[c,lags] = xcorr(genCode1,'normalized');
plot(lags,c)
%% 基带信号
basebandSignal1=zeros(1,7500);
tau=1; %时延/码片
alpha=0.5;%衰落系数
basebandSignal0=basebandSignal;%直达信号
basebandSignal1(1+tau:end)=basebandSignal(1:end-tau);
basebandSignal1(1:tau)=0;
basebandSignal_r=basebandSignal0+alpha*basebandSignal1;
[c,lags] = xcorr(basebandSignal_r,basebandSignal0,'normalized');
plot(lags,c)
set(gca,'XLim',[-5 5]);%X轴的数据显示范围
%% 接收信号
[c,lags] = xcorr(basebandSignal, promptCode,'normalized');
plot(lags,c)
set(gca,'XLim',[-30 30]);%X轴的数据显示范围
%%
% 设置参数
clear
fs = 10000;                   % 采样率
fc = 100;                    % 信号载频
fd = 5;                      % 多普勒频移
doppler_bw = 0;              % 多普勒带宽
num_paths = 2;               % 多径数量
path_delays = [0, 0.5] * 1e-6;% 多径时延
path_gains = [0,-10];          % 多径增益
t = 0:1/fs:1;                % 时间向量
M=4;                         % 调制指数
navmsgFreq=50;               % 码频率

% 生成原始信号
bit_signal = round(rand(1,1000));                                   % 输入数据流
t = 0:1 /fs:1;%时间向量
phaseCache = navmsgFreq * t;                                        % 导航电文相位
genNavmsg = bit_signal(floor(phaseCache(1:(end-1))) + 1);           % 信号导航电文
figure
plot(t(1:end-1),genNavmsg)
phaseCache = 2 * pi * fc * t;                                       % 载波相位
genCarr = exp(1i * phaseCache(1:(end-1)));                          % 信号载波
x= genNavmsg.*genCarr;
figure
plot(t(1:end-1),x);
N = length(x);
f = (-N/2:N/2-1)*(fs/N);
Y = fftshift(fft(x));
figure
subplot(3,1,1);
plot(f, abs(Y));
title('发射信号频谱');
% 生成多径信道
max_path_delay = max(path_delays);
num_samples = length(t);
channel = zeros(num_paths, num_samples);
for i = 1:num_paths
    delay_samples = round(path_delays(i)*fs);
    channel(i, delay_samples+1:end) = 10^(path_gains(i)/20) ...
        * exp(-1i*2*pi*doppler_bw*(0:num_samples-delay_samples-1)/fs);
end
h = sum(channel);
FadedSignal = filter(h,1,x);
% 多普勒频扩
y_end = conv(x, h);
figure
plot(t(1:end-1),FadedSignal);

bit_moded3 = qamdemod(y_end, M, 'OutputType', 'bit'); % 解码
% 频谱分析
N = length(y_end);
f = (-N/2:N/2-1)*(fs/N);
Y = fftshift(fft(y_end));
figure
subplot(3,1,2);
plot(f, abs(Y));
title('接收信号频谱');
% 绘制信道频率响应
H = fftshift(fft(h, N));
subplot(3,1,3);
plot(f, abs(H));
title('信道频率响应');
basebandSignal = FadedSignal .* genCarr;
D=sqrt(real(basebandSignal).^2+imag(basebandSignal).^2);
%%
for i=1:length(D)
    if D(i)>1
        D(i)=1;
    end
end
figure
plot(t(1:end-1),D,'.')
%% 观测载波信号
genSignal = genCarr.* genNavmsg;
t=(0:genSampleCount) / sampleFreq;
t=t(1:1000);
y_end=genSignal(1:1000);
plot(t,y_end);
%% 可视化
hold on
plot(7/10,mean(D),'*');
set(gca,'XLim',[0 1]);%X轴的数据显示范围
set(gca,'YLim',[0 0.3]);%Y轴的数据显示范围
%% 积分测试
for m=1:1000
    t=0:1/m:pi;
    y_end=sin(t);
    d(m)=mean(y_end);
end
plot(d)
%%
    delta_tau0=0:0.00001:2;
    fc=1575.42e6;%载波频率
    L=1;%时延变化范围（码片）
    %delta_t=-L:L/2000:L;%抖动时延
    %delta_tau1=delta_tau0+delta_t;
    delta_tau1=delta_tau0;
    alapha1=0.5;%衰落系数
    delta_fai_1=delta_tau1*2*pi*fc*1e-6;%多径相位时延变化
    fai1=mod(delta_fai_1,2*pi);
    fai_c=atan(alapha1*sin(delta_fai_1)./(1+alapha1*cos(delta_fai_1)));%载波跟踪环路相位
    fai=delta_fai_1+fai_c;%相位
    fai00=mod(fai,2*pi);
    Tc=1;%C/A码片宽度
    d=Tc;%早迟间隔
    tau_L=(1+alapha1*cos(fai))/2*d;
    tau_H=(alapha1*d*cos(fai))/2+Tc-d/2;
    error1=(delta_tau1.*alapha1.*cos(fai)./(1+alapha1.*cos(fai)));
    figure
    plot(delta_tau1,error1,'*')
    figure
    plot(delta_tau1,fai1,'-')
    hold on
    plot(delta_tau1,fai_c,'-')
    plot(delta_tau1,fai00,'-')
    plot(delta_tau1,error1,'-')
    xlim([0.3,0.302]);
    xlabel('时延 / 码片');
    ylabel('码多径误差/码片');
    %% 算积分
% 选择积分指令int()，具体用法可查阅帮助文档
% 确定积分元x,将最后的积分结果保存在f中
clear
syms f x a;
f = int(a*cos(x)./(2-a*cos(x)),x);  % 计算不定积分

pretty(f)  % 格式化显示积分结果
%%
clear
a=0.5;
x=pi+0.00001:0.1:3*pi;
y_end=-x-(2*log((a - 2*cos(x) + sin(x)*(a^2 - 4)^(1/2))./(a*cos(x) - 2)))/(a^2 - 4)^(1/2);
q1=y_end(end)-y_end(1)
a=0.5;
x=pi+0.00000000001:0.01:3*pi-0.000000000001;
y_end=x - log((a + cos(x) + sin(x).*(a^2 - 1)^(1/2))./(a.*cos(x) + 1))./(a^2 - 1)^(1/2);
x_begin=x(1);
x_end=x(end);
y_1=x_begin-log((a+cos(x_begin)+sin(x_begin)*sqrt(a^2-1))./(a*cos(x_begin)+1))./(sqrt(a^2-1));
y_end=x_end-log((a+cos(x(end))+sin(x_end)*sqrt(a^2-1))./(a*cos(x_end)+1))./(sqrt(a^2-1));
q2=y_end(end)-y_end(1)

q2/q1
%% 理论结果
clear
a=0.5;
x_en=0.001:0.5:16*pi;
for ll=1:length(x_en)
    x=10*pi+0.00000000000001:0.001:10*pi+x_en(ll);
    y2=a*cos(x)./(1+a*cos(x));
    x_begin=mod(x(1)-pi,2*pi)+pi;
    x_end=mod(x(end)-pi,2*pi)+pi;
    aa=x(1)-x_begin;
    bb=-x(end)+x_end;
    d=(aa+bb)/2/pi;
    y_1=log((a+cos(x_begin)+sin(x_begin)*sqrt(a^2-1))./(a*cos(x_begin)+1))./(sqrt(a^2-1));
    y_end=log((a+cos(x(end))+sin(x_end)*sqrt(a^2-1))./(a*cos(x_end)+1))./(sqrt(a^2-1));
    y_mean=(-y_end+y_1+x_end-x_begin+d*0.9688)./(x(end)-x(1));
    erro1(ll)=mean(y2);
    erro2(ll)=y_mean;
end
plot(x_en/pi,abs(erro2),'r-');
hold on
plot(x_en/pi,abs(erro1),'>')
erro1(end)
%%
plot(x/pi,y_yhs);
hold on
plot(x/pi,y2*3);
figure
plot(x_en/pi,abs(error));
%% 理论
clear
a=0.5;
    x=pi+0.00000000000001:0.001:3*pi-0.00000000000001;
    y2=a*cos(x)./(1+a*cos(x));
    x_begin=x(1);
    x_end=x(end);
    y_yhs=x-log((a+cos(x)+sin(x)*sqrt(a^2-1))./(a*cos(x)+1))./(sqrt(a^2-1));
    y_yhs_1=log((a+cos(x_begin)+sin(x_begin)*sqrt(a^2-1))./(a*cos(x_begin)+1))./(sqrt(a^2-1));
    y_yhs_end=log((a+cos(x(end))+sin(x_end)*sqrt(a^2-1))./(a*cos(x_end)+1))./(sqrt(a^2-1));
    aa=x(1)-x_begin;
    bb=-x(end)+x_end;
    d=-(aa+bb)/2/pi;
    y_1=x_begin-log((a+cos(x_begin)+sin(x_begin)*sqrt(a^2-1))./(a*cos(x_begin)+1))./(sqrt(a^2-1));
    y_end=x_end-log((a+cos(x(end))+sin(x_end)*sqrt(a^2-1))./(a*cos(x_end)+1))./(sqrt(a^2-1));
    y_mean=(y_end-y_1)./(x(end)-x(1))
plot(x/pi,y_yhs);
%plot(x/pi,y2*3);
%%
clear
alapha1=0.5;
t=0:0.0001:2;
fc=1525;
fai=2*pi*fc.*t;
y_end=t.*alapha1.*cos(fai)./(1+alapha1.*cos(fai));
plot(t,y_end)